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ABSTRACT 

Graham et al. 2015 discovered a supermassive black hole binary (SMBHB) candidate 
and identified the detected 5.2 yr period of the optical variability as the orbital period 
of the binary. Hydrodynamical simulations predict multiple periodic components for 
the variability of SMBHBs, thus raising the possibility that the true period of the 
binary is different from 5.2 yr. We analyse the periodogram of PG1302 and find no 
compelling evidence for additional peaks. We derive upper limits on any additional 
periodic modulations in the available data, by modeling the light-curve as the sum of 
red noise and the known 5.2 yr periodic component, and injecting additional sinusoidal 
signals. We find that, with the current data, we would be able to detect with high 
significance (false alarm probability < 1%) secondary periodic terms with periods 
in the range predicted by the simulations, if the amplitude of the variability was 
at least ~0.07mag (compared to 0.14mag for the main peak). A three-year follow¬ 
up monitoring campaign with weekly observations can increase the sensitivity for 
detecting secondary peaks «50%, and would allow a more robust test of predictions 
from hydrodynamical simulations. 


1 INTRODUCTION 

It is well established that massive galaxies harbour super¬ 
massive black holes in their centres, with masses tightly cor¬ 
related with global properties of the host galaxy (Kormendy 
& Ho 2013). According to cosmological models of hierar¬ 
chical structure formation, galaxies grow through frequent 
mergers (Springel, Di Matteo & Hernquist 2005; Robertson 
et al. 2006). The formation of supermassive black hole bi¬ 
naries (SMBHB) is a natural consequence of the above pro¬ 
cess. Combining the high expected rates of galaxy mergers 
with the expectation that SMBHBs spend a large fraction of 
their lifetimes at close separation (Haehnelt & Kauffmann 
2002), SMBHBs at sub-parsec separations should be fairly 
common, despite the scarcity of observational evidence. 

Recently, Graham et al. (2015, hereafter G15) reported 
the detection of strong periodic variability in the optical flux 
of quasar PG1302-102. PG1302 is a bright (median V-band 
magnitude ~15), radio-loud quasar at redshift z = 0.2784, 
with inferred black hole (BH) mass of 1O®'^“® "'M 0 . The 
light curve in optical bands shows a quasi-sinusoidal modu¬ 
lation, with a best-fit period of (5.2 ± 0.2) yr and amplitude 
of «0.14mag. G15 suggest that PG1302 may be a SMBHB 
at close separation (~ 0 . 01 pc), interpreting the observed pe¬ 
riodicity as the (redshifted) orbital period of the binary. 

Theoretical work on circumbinary disks predicts that 
SMBHBs can excite periodic enhancements of the mass ac¬ 
cretion rate that could translate into periodic luminosity 
enhancements, not only at the orbital period, but also on 
longer and shorter timescales. More specifically, if the binary 


is embedded in a thin disk, the gas will be expelled from the 
central region due to torques exerted by the binary, creating 
a cavity (Artymowicz & Lubow 1994). Several hydrodynam¬ 
ical simulations (Hayasaki, Mineshige & Sudou 2007; Mac- 
Fadyen & Milosavljevic 2008; Noble et al. 2012; Shi et al. 
2012; Roedig et al. 2012; D’Orazio, Haiman & MacFadyen 
2013; Farris et al. 2014) indicate that the interaction of the 
individual BHs with the inner edge of the accretion disk can 
pull gaseous streams into the cavity, resulting in periodic 
modulation of the accretion rate at timescales correspond¬ 
ing to «1/2 and 1 times the binary’s orbital period. For high 
BH mass ratios {q = M 1 /M 2 > 0.3), the cavity is lopsided, 
leading to the formation of a hotspot in the accretion disk. 
The strongest modulation in the accretion rate is observed 
at the orbital period of the overdense region, a few (~3-8) 
times the orbital period of the binary. 

These results imply that the observed 5.2 yr period in 
PG1302 may not reflect the orbital period of the binary. 
If the orbital period of the binary is shorter/longer than 
the period of optical valiability (hereafter, optical period, 
topt), there are major implications for the expected quasar¬ 
binary fraction, and the probability of detecting SMBHBs 
(Haiman, Kocsis & Menou 2009), as well as the physics of 
the orbital decay. D’Orazio et al. (2015) showed that, under 
the assumption that the optical period corresponds to the 
longer period of the hotspot in the accretion disk, PG1302 
would be in the gravitational inspiral regime. This would 
confirm that SMBHBs can produce bright electromagnetic 
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emission even at these late stages of the merger (Noble et al. 
2012; Farris et al. 2015). 

The significance of the above considerations motivates 
us to search for multiple periodicity in the optical variabil¬ 
ity of PG1302. In this paper, we search for secondary peaks 
in the Lomb-Scargle periodogram. We assess the statistical 
signihcance of the identified peaks by generating mock time 
series that show correlated noise, as expected for quasars 
(Kelly, Bechtold & Siemiginowska 2009). The false alarm 
probability of the most significant secondary peak identified 
in our analysis is 6%, below a reliable detection threshold. 
We set upper limits on the amplitude of secondary sinusoid 
variations that could be detectable with high significance, 
in the presence of correlated noise, as well as the already 
known 5.2 yr periodic signal. Although the current data do 
not allow the detection of weak secondary periodic terms, 
the addition of three years of observations, can improve the 
sensitivity up to a factor of 2, lowering the detection thresh¬ 
old to ~0.04mag. 


2 DATA ANALYSIS 

We investigate the possibility of multiple periodic terms in 
the photometric variability of PG1302 by analysing the light 
curve published in G15^. The light curve consists of data 
from: (1) the Gatalina Real-time Transient Survey (tele¬ 
scopes GSS and MLS; Drake et al. 2009; Mahabal et al. 
2011); (2) the Lincoln Near-Earth Asteroid Research (LIN¬ 
EAR; Sesar et al. 2011) program; (3) the All Sky Automated 
Surveys (ASAS; Pojmanski 1997) and; (4) other archival 
data (Garcia et al. 1999; Eggers, Shaffer & Weistrop 2000), 
providing a baseline of ~20yr. The different datasets were 
calibrated to account for the differences between the various 
photometric systems, as detailed in G15. 

For our analysis, we use the generalised version of the 
Lomb-Scargle (LS) periodogram (Lomb 1976; Scargle 1982; 
Zechmeister & Kiirster 2009), a standard method for de¬ 
tecting periodicity in unevenly sampled time series^. We 
compute the periodogram for 1000 trial frequencies on a 
uniform logarithmic grid, spanning from fmin = ijTdata to 
fmax = 1/(2AT), where Tdata is the baseline of the light 
curve and AT is the median time interval between succes¬ 
sive data points. 

For each peak in the periodogram, we then calculate the 
probability that a peak of equal power can arise by chance 
(false alarm probability; hereafter FAP). For this purpose, 
we generate mock time series that mimic the optical variabil¬ 
ity of PG1302. Quasars show correlated stochastic variabil¬ 
ity, best described as a damped random walk (Kelly, Bech¬ 
told & Siemiginowska 2009). At high frequencies, the power 
spectral density decreases with frequency {PSD ~ 1//^, 
or red noise), whereas, at low frequencies, it becomes flat 
{PSD ~ /°, or white noise). Furthermore, PG1302 ex¬ 
hibits strong sinusoidal variability with period of 1, 884±88 d 
(intentified by wavelet and autocorrelation function (ACF) 


^ We are grateful to M. Graham for providing us with their cal¬ 
ibrated photometric data in electronic form. 

^ Throughout the analysis, we make use of the astroML python 
package (Vanderplas et al. 2012; Ivezic et al. 2014). 


based techniques as detailed in G15) and amplitude of 0.14 
mag. 

We first identify the best-fit model for PG1302 as fol¬ 
lows. We generate evenly sampled stochastic time series with 
a dense temporal resolution {dt = 2h) and peak-to-peak am¬ 
plitude A, and with a power spectrum ~ 1//^ (Timmer & 
Koenig 1995). We add a sine wave, with properties as in 
G15, and subsample to match the observation times in the 
light curve of PG1302. We calculate the average LS peri¬ 
odogram from 50 realisations. Next, we subtract this peri¬ 
odogram from the observed LS periodogram of PG1302 and 
calculate the residuals. We repeat the above process, vary¬ 
ing the amplitude A between 0.5 to 2 times the peak-to-peak 
amplitude of the light curve {rrimax — mmin) of PG1302, and 
determine the amplitude A that minimises the residuals. 

With the best-fit amplitude A, we generate 10,000 re¬ 
alisations of red noise with periodic variability, as before, 
and calculate the LS periodogram for each realisation. At 
each trial frequency, we calculate the FAP by comparing 
the power in the PG1302’s periodogram to the distribution 
of power at the same frequency in the periodograms of the 
mock time series. We dehne FAP/, the false alarm proba¬ 
bility per frequency as the fraction of the 10,000 simulated 
periodograms with power exceeding the value observed for 
PG1302. 

Based on the results from hydrodynamical simulations 
discussed above, we search for secondary peaks around the 
main 5.2yr peak, i.e. between (/min,/max) = (2.37, 3.8) nHz 
and (9.5, 92) nHz®. The frequency ranges together contain a 
total of Ntot = 526 frequency bins. Since these bins have 
non-trivial correlations, we compute the effective number 
of independent bins, Neg, for each FAP/, by: (i) applying 
the same analysis as above, but replacing the PG1302 data 
with each individual simulated mock periodogram, and (ii) 
determining the fraction of the 10,000 mock periodograms 
that show FAP / (or lower) at any frequency between /min < 
/ < /max. For FAP/ = 10“^, we find N^g = 80, and the total 
FAP= A'effXFAP/=0.008.'‘ 


3 RESULTS 

If the optical variability of PG1302 is the superposition of 
multiple periodic terms, it is expected that significant sec¬ 
ondary peaks will appear in the periodogram. Moreover, fol¬ 
lowing the predictions of hydrodynamical simulations, we 
expect the peaks roughly at ratios 1:2 (1 and 1/2 times the 
orbital period) or from 1:3:6 to 1:8:16 (hotspot period, or¬ 
bital period and half of the orbital period of the binary). 

More specifically, with the limitations imposed by the 
currently available data, we were able to explore three pos¬ 
sibilities: (A) If the optical period tapt is the binary orbital 
period tbin, is there a secondary peak near 0.5tbin? (B) If 
topt is O.Stbin, is there a periodic component at tbin? (C) If 
topt is the long orbital period « (3 — 8)tbin of the hotspot in 
the disk, as hypothesised by D’Orazio et al. (2015), is there 
a secondary period corresponding to « O.Stbin or tbin? Note 
that in scenario (A), periodicity is predicted at ~ (3 —8)tbin- 

® These ranges are illustrated by the shaded areas in Fig. 3 
* Nog is decreased to ^ 60 and ss 30 for FAP/ = 10“® and 
FAP/ = 10“®, respectively. 
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Figure 1. LS periodogram of PG1302 (blue) and of a sine wave 
sampled at the observation times of PG1302 (red). The embedded 
plots illustrate the respective light curves, spanning a baseline of 
~20 yr. The irregular sampling of the data produces the artificial 
secondary peaks seen in the red curve. 

However, the 20-yr baseline of the light curve does not allow 
us to probe this low-frequency range. 

The detection of weaker periodic terms in the presence 
of a known significant periodic component is in general chal¬ 
lenging. The periodogram of an unevenly sampled time se¬ 
ries is the convolution of the actual power spectrum of the 
signal and the window function, defined by the irregular time 
sampling. A periodic signal can introduce aliasing peaks, 
and hide the secondary peaks. To illustrate this, in Fig. 1, we 
show the LS periodogram of PG1302, superimposed with the 
periodogram of a pure sine wave, with period of 5.2 yr, sam¬ 
pled at the observation times of PG1302. The periodogram 
of the sine wave shows a set of artificial secondary aliasing 
peaks, which coincide with many of the observed features in 
the periodogram of PG1302.® 

We analysed the LS periodogram of PG1302 searching 
for secondary peaks, focusing on the frequencies of interest, 
within a factor of ~10 around the peak of the main pe¬ 
riod. We assessed the statistical significance of the observed 
secondary peaks using the method discussed in § 2 above. 
In Fig. 2, we present the periodogram of PG1302 with the 
average and the maximum periodogram power from 10,000 
realisations of the simulated variability, with the latter cor¬ 
responding to our FAP=1% significance threshold. In the 
bottom panel, we show both FAP/ and the total FAP at 
each frequency. We identify a secondary peak at ~40 nHz 
(~300d), with the lowest FAP (< 0.008). This peak is most 
likely artificial, since it coincides with one of the strongest 
aliasing peaks demonstrated in Fig. 1. The most significant 
peak that does not overlap with aliasing peaks appears at 
~77.5nHz (~155d), with a FAP=0.06 (FAP/=0.001). 

We next derive upper limits for putative secondary pe¬ 
riodic terms with the current data. For this purpose, we 
simulate time series as before, but with an additional sine 
wave component. We identified the minimum amplitude of 
the second component at which the power in the new peak 

® The aliasing peaks are predictable, and could be used to better 
determine the amplitude and frequency of the 5.2 yr peak. 



Frequency[Hz] 


Figure 2. Top panel: LS periodogram of PG1302 (blue), super¬ 
imposed with the average of 10,000 realisations of a model that 
includes stochastic red noise and a 5.2 yr-sinusoid. The maximum 
power among these 10,000 realisations is also shown (black), cor¬ 
responding to 1% false alarm probability threshold. Bottom panel: 
False alarm probability per frequency (red) and total FAP (black) 
of PG1302’s periodogram, as a function of frequency. Note that 
the y-axis in the bottom panel decreases upwards. 

in the LS periodogram exceeded the FAP < 1% detection 
threshold at least 90% of the time. We repeated this process 
for different frequencies within the factor of 10 range of in¬ 
terest. In Fig. 3, we present the minimum relative amplitude 
(corresponding also to V-band magnitude) of a secondary si¬ 
nusoid that would be detectable at different frequencies. The 
shaded areas in this figure indicate frequencies relevant to 
the three possibilities (A)-(C) discussed above. 

As Fig. 3 shows, we can set weak limits for scenario (A); 
a second periodic component would need a higher amplitude 
than the one identified in G15 to be detectable. The reduced 
sensitivity in this frequency range is reasonable, since the 
light curve is well sampled only in the last ~9yr, hindering 
the detection of a ~ 10 yr period. Moreover, at this frequency 
range the effect of the red noise is significant. For scenario 
(B) and for the majority of frequencies relevant to scenario 
(G), the secondary term needs to have amplitnde compara¬ 
ble to the main 5.2 yr modulation to be detectable. Finally, 
for a few specific frequencies relevant to scenario (G), we can 
probe secondary sinnsoids with amplitudes down to 25% of 
the amplitude of the main sinusoid (~ 0.04mag). For the 
peak with the lowest FAP (155 d), the minimum detectable 
amplitude would be 50% of the main sinusoid (0.07mag). 


4 DISCUSSION 

4.1 Quasar variability models 

We have thus far modelled the variability of PG1302 us¬ 
ing a fixed power spectrum 1//^. Since quasar variability is 
described by a damped random walk (i.e. with the power 
spectrum flattening at low frequencies), white (or pink) 
noise may describe the quasar variability more accurately 
at low frequencies. To assess how this impact our conclu¬ 
sions, we repeated the analysis described above but with 
stochastic noise exhibiting 1/f (pink noise) and flat (white 
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Figure 3. The minimum amplitude of a secondary sinusoid term 
that would be detectable at each frequency. The y axis shows the 
amplitude relative to that of the main 5.2 yr peak (left), and also 
in V mag (right). The three curves correspond to noise models 
with three different power-law power spectra: oc 1//^ (brown), 
oc 1// (olive), and constant white noise (dark blue). Frequency 
ranges relevant to different scenarios are highlighted with light 
pink (A), light green (B) and light yellow (C). Vertical lines mark 
specific frequencies of interest (see text). 

noise) power spectra. Our conclusions remain unchanged for 
these different noise models. More specifically, we identify 
the same peak at ~40nH (~300d), which we discard as 
aliasing with FAP< 1% both for white and pink noise; we 
also find the peak at ~77.5nH (~155d) with FAP=0.2 and 
0.1 for pink and white noise, respectively. The upper limits 
for the different noise models are also shown in Fig. 3. For 
white noise, the minimum detectable amplitude is almost 
constant (~0.08mag) for all frequencies, but with fluctua¬ 
tions caused by the irregular sampling similarly to the red 
noise case. The results are again similar for pink noise, with 
the sensitivity falling in-between the red and white noise 
cases. 

4.2 Historic Data 

Historic photometry of PG1302 is available from the digiti¬ 
sation of old photographic plates, as part of the Digital Ac¬ 
cess to a Sky Century @ Harvard (DASCH; Grindlay et al. 
2009) project. The measurements were in unfiltered B pho¬ 
tographic magnitude, which is close to wide-band Johnson 
B magnitude and were calibrated with the AAVSO APASS 
catalog (Henden et al. 2012). We extracted the B-band light 
curve from the online database® and excluded measurements 
with high astrometric errors or with magnitudes within 0.5 
of the limiting magnitude, and data points that were flagged 
as plate defects. We analysed the periodogram, as before, 
assuming red noise variability (without the periodic compo¬ 
nent) and finding the best-fit model. 

In Fig. 4, we present the B-band light curve and the LS 
periodogram, along with the average and maximum power of 
the LS periodogram from 10,000 realisations of red noise. All 
of the identified peaks (including the 5.2 yr period from G15) 


Figure 4. Top Panel: B-band light curve of PG1302 from the 
DASCH project. Bottom Panel: LS periodogram of the B-band 
light curve (blue), superimposed with the average best fit 1//^ 
noise model (red) and the maximum power of the periodogram 
from 10,000 realisations of the noise (black), which corresponds 
to our significance threshold. 

are below the significance threshold. This result, however, 
must be viewed in light of the photometric accuracy of the 
DASGH dataset. The average photometric error is 0.18 mag, 
exceeding the 0.14 mag amplitude of the periodic variability 
in G15. Interestingly, the two most significant peaks, within 
the interval of interest, are observed at similar frequencies 
(~42 nH and ~76 nH) as the peaks in the periodogram from 
the G15 data, although their significance remains low. 

4.3 Future Data 

We investigated how the sensitivity for detection of sec¬ 
ondary peaks would increase with three years of weekly 
follow-up observations. In order for PG1302 to be observ¬ 
able throughout the year, a network of telescopes is required 
(e.g., Las Gumbres Observatory^). Given that PG1302 is 
a bright quasar, even small telescopes can observe it with 
relatively small photometric errors. To allow for a flexible, 
realistic (affected by weather conditions, resources etc.) ob¬ 
servation schedule as well as to avoid strong aliasing from 
an exactly periodic time sampling, we generate a hypothet¬ 
ical follow-up data set, consisting of one randomly chosen 
day each week, and allowing an additional scatter of ±3 h in 
the observation times. We also assume photometric errors of 
lOmmag, comparable to the ones in G15. 

We analyse the extended light curve as before: we gen¬ 
erate 10,000 time series that exhibit red noise variability 
with a periodic component as in G15, sampled at the obser¬ 
vation times (including both the existing and future data) 
and calculate the LS periodogram. We set our significance 
threshold to the maximum power of the periodogram from 
the 10,000 realisations (FAP<1%). Next, as done above, we 
include a secondary sinusoidal term and identify the min¬ 
imum amplitude that would result in a peak with power 
above the threshold 90% of the time. 


http://dasch.rc.fas.harvard.edu/lightcurve.php 


^ http://lcogt.net/ 
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Figure 5. Minimum relative amplitude (also in V mag) of a 
secondary sinusoid versus frequency of secondary sinusoid that 
could be detectable in the current data (brown), in data extended 
with three years of weekly follow-up observations (olive) and in 
three years of future data only (dark blue). The shaded regions 
show the frequencies of interest, with color coding as before. 


Fig. 5 shows the upper limits for detection with the ex¬ 
tended data, compared to the upper limits with only the 
current data. For low frequencies (e.g., scenarios (A) and 

(B) ), the sensitivity does not improve with the addition of 
three-years of data. A longer monitoring campaign is re¬ 
quired for this purpose. For frequencies relevant to scenario 

(C) , the improvement varies from marginal to a factor of 
2, depending on the precise frequency. Furthermore, we cal¬ 
culate upper limits based only on the future data and find 
that the sensitivity for detection of a secondary term can 
reach as low as 20% of the amplitude of the main compo¬ 
nent (~0.03 mag) at some frequencies. This suggests that a 
multi-year follow-up with frequent sampling, by itself, could 
significantly improve the sensitivity to detect high-frequency 
(~ 1 yr) peaks; 


4.4 Implications for SMBH Binary Models 

The lack of a significant secondary peak can place con¬ 
straints on the physical parameters of the SMBHB PG1302. 
In particular, models that predict two peaks, comparable 
in height, and separated in frequency by a factor of few, 
can be ruled out. While this is not a generic prediction of 
hydrodynamical simulations, for specific mass ratios within 
the range of 0.25 < q < 0.45, two near-equal peaks are pre¬ 
dicted in the periodograms of the mass accretion rate onto 
the BHs Farris et al. (2014, see their Fig. 9). The specific 
value of the disfavoured ratio q depends on disk parame¬ 
ters (viscosity, temperature) - and also on how one defines 
the accretion rate. If the accretion rate is defined based on 
the total mass accreting into the central cavity, this q value 
would be somewhat above q = 0.5 (see Fig. 6 in D’Orazio, 
Haiman & MacFadyen 2013). 


5 CONCLUSIONS 

The presence of multiple periodic components in the vari¬ 
ability of SMBHBs is predicted by hydrodynamical simula¬ 


tions of circumbinary disks. The detection of multiple pe¬ 
riodicity would provide additional indication that PG1302 
is a SMBHB. We analysed the LS periodogram of the op¬ 
tical light-curve of the SMBHB candidate quasar PG1302 
to search for multiple periodic components, and assessed 
the statistical signihcance of secondary peaks by simulat¬ 
ing the variability of PG1302. Our analysis returned a peak 
with FAP=6%, below our detection threshold. By injecting 
fake secondary sinusoids into our models, we found that the 
current data would only allow the detection of secondary 
periodic modulations comparable in amplitude to the main 
peak. Future observations (for a few years, at weekly sam¬ 
pling), and/or a more sophisticated analysis, which can mit¬ 
igate the effects of aliasing, could signihcantly improve the 
detection capability and help uncover secondary peaks ex¬ 
pected to be produced by SMBH binaries. 
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